Flow velocimeter system

ABSTRACT

While a light sheet is generated at a designated region, images of fluid flowing through the designated region are formed at different times. For an inspection region of the plurality of inspection regions defined in the images that has a degree of difference exceeding a threshold between the local flow velocity vector v(a, b, T) at a certain time T and a reference flow velocity vector v(a, b, T ± ) at times T ±  that are different from the certain time T, the flow velocity vector v(a, b, T) at the reference time T is corrected with the reference flow velocity vector v(a, b, T ± ).

BACKGROUND OF THE INVENTION

Field of the Invention

The present invention relates to flow velocimeter systems.

Description of the Related Art

Known flow velocimetry includes a particle image velocimeter(hereinafter called a PIV) that is configured to irradiate smallparticles mixed in fluid with light such as laser light to acquire animage of scattered light therefrom, thus finding a moving distance of aparticle or a group of particles from images of particles acquired atslightly different times and so measuring the velocity of fluid (flowvelocity) in the two-dimensional plane.

In the velocity field of the PIV, vectors (noise) having a very largeerror are present. Such noise is generated because particles flow outfrom or into a light sheet due to the velocity component perpendicularto the light sheet or the voracity component, thus failing to find acorrelation between the images. Such generation of the noise cannot beavoided, meaning that the velocity obtained includes a very large error.

To solve such a problem, Japanese Patent Application Laid-Open No.2004-20385 (Patent Document 1) proposes a technique of performingFourier transformation of a time-series velocity at each velocitydefining point in flow velocity vector distributions acquired from twoimages that are temporally continuous, to which cutoff filtering orlow-pass filtering is applied in the frequency space for inverse Fouriertransformation, thus obtaining a final time-series velocity of thefluid.

SUMMARY OF THE INVENTION

Such a conventional technique, however, is configured to performprocessing such as Fourier transformation to a time-series velocity ateach velocity defining point as one continuous function, and so has tobe modified better from the viewpoint of improving the accuracy toacquire the time-series velocity of the fluid, because while an error ata noise part having a large error can be decreased, but an error at aflow velocity vector part having a small error increases unfortunately.

Then the present invention aims to provide a system for PIV that iscapable of improving the measurement accuracy of time-series or changingmode of the fluid velocity.

A flow velocimeter system of the present invention comprises: alight-sheet generation device; an imaging device; and an image analysisdevice including a processor, the image analysis device being configuredto measure, based on images of fluid flowing through a designated regionthat are formed a plurality of times at different times by the imagingdevice while the light-sheet generation device generates a light sheetcontinuously or intermittently at the designated region, a time-serieschanging mode of a local flow velocity vector of the fluid at each of aplurality of inspection regions defined at the images in accordance witha cross correlation method. The image analysis device is configured to,for an inspection region of the plurality of inspection regions that hasa degree of difference exceeding a threshold between the flow velocityvector at a certain time and a reference flow velocity vector as theflow velocity vectors at one or a plurality of times that is differentfrom the certain time, correct the flow velocity vector at the certaintime with the reference flow velocity vector.

According to the flow velocimeter system of the present invention, alight sheet is generated at the designated region, and images of fluidflowing through the designated region are formed at a plurality ofdifferent times. For each of a plurality of inspection regions making upthe images, when a flow velocity vector at a certain time is differentfrom a flow velocity vector (reference flow velocity vectors) at one ora plurality of times that is different from the certain time beyond athreshold, the flow velocity vector at the certain time is corrected. Asa result, a flow velocity vector that is likely to be noise differentfrom the original form can be corrected closer to the original form.Then, the accuracy of measurement of time-series or changing mode of thefluid velocity can be improved.

In the flow velocimeter system in one aspect of the present invention,the image analysis device is configured to correct brightness of allpixels making up a predetermined range that is at least a part of theimages based on brightness of a pixel as a part making up thepredetermined range, and then calculate a flow velocity vector at eachof the plurality of inspection regions.

According to the thus configured flow velocimeter system, brightness ofall pixels making up a predetermined range that is at least a part ofthe formed image is corrected. This enables clear discrimination of thebrightness resulting from the border between particles included in thefluid or making up the fluid and the background light and the brightnessresulting from the scattered light of the particles. That is, themeasurement accuracy of a time-series correlation of the regionincluding the particles in the formed image can be improved, andaccordingly the time-series measurement accuracy of the fluid velocitycan be improved.

In the flow velocimeter system of the present invention, the imageanalysis device is configured to execute fast Fourier transformationprocessing at each of the plurality of inspection regions, thuscalculating the flow velocity vectors, and when there is an overlappedregion in one inspection region and in another inspection region amongthe plurality of inspection regions, to apply a result of the fastFourier transformation processing at the other inspection region to theoverlapped region for the fast Fourier transformation at the oneinspection region.

According to the thus configured flow velocimeter system, thecalculation resource can be saved for the fast Fourier transformationprocessing (FFT processing) in the overlapped region between theplurality of inspection regions. This means that the fluid velocityhaving a high pixel density of the image also can be measuredeffectively in a time-series manner, and so the measurement accuracy canbe further improved.

In the flow velocimeter system of the present invention, the imageanalysis device uses a value that is determined in accordance with thereference flow velocity vector as the threshold.

According to the thus configured flow velocimeter system, the thresholdcan be determined adaptively corresponding to the reference flowvelocity vector, and so the time-series measurement accuracy of thefluid velocity can be further improved.

In the flow velocimeter system of the present invention, the imageanalysis device is configured to, when occupancy of different flowvelocity vectors among the reference flow velocity vectors at theplurality of times exceeds a predetermined value, increase the number ofthe reference flow velocity vectors so that the occupancy becomes thepredetermined value or less. The different flow velocity vectors arereference flow velocity vectors at one time of the plurality of times,which have a degree of difference from reference flow velocity vectorsat another time that exceeds the threshold.

According to the thus configured flow velocimeter system, the number ofreference flow velocity vectors (imaging times) are increased as neededso as to reduce the influences from reference flow velocity vectors thatare likely to correspond to noise. Then, the flow velocity vector ateach time can be corrected suitably based on the reference flow velocityvector that is likely not to correspond to noise or to correspond to thetrue value, and so the time-series measurement accuracy of the fluidvelocity can be improved.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 schematically illustrates the system as a whole.

FIG. 2 illustrates the relationship between the exposure duration of animaging device and timing of laser irradiation.

FIG. 3 is a flowchart to acquire a time-series velocity of fluid.

FIG. 4 illustrates the relationship between the image size of a particleand an error of cross correlation.

FIG. 5 describes particle image correction processing.

FIG. 6A illustrates the range where flow velocity vectors can beacquired when the particle image correction processing is not performed.

FIG. 6B illustrates the range where flow velocity vectors can beacquired when the particle image correction processing is performed.

FIG. 7A describes a first inspection region located at the upper leftcorner.

FIG. 7B describes the first inspection region that is displaced by onein the x-axis direction.

FIG. 7C describes the first inspection region that is displaced by onein the y-axis direction.

FIG. 7D describes a typical first inspection region.

FIG. 8A describes a second inspection region.

FIG. 8B describes the calculation of a correlation function.

FIG. 8C describes the calculation of a flow velocity vector.

FIG. 9 describes the brightness distribution at each inspection regionand FFT processing to be performed for the brightness distribution.

FIG. 10A, FIG. 10B, and FIG. 10C describe a part where FFT processingcan be omitted for the brightness distribution at the inspectionregions.

FIG. 11 is a flowchart of time-series error correction processing.

FIG. 12A describes the correction of each flow velocity vector in thetime-series error correction processing.

FIG. 12B describes a change in window length in the time-series errorcorrection processing.

FIG. 13 describes window setting in the time-series error correctionprocessing.

FIG. 14A illustrates a time-series changing mode of original data and aflow velocity vector of fluid before the time-series error correctionprocessing.

FIG. 14B illustrates a time-series changing mode of original data and aflow velocity vector of fluid after the conventional time-series errorcorrection processing.

FIG. 14C illustrates a time-series changing mode of original data and aflow velocity vector of fluid after the time-series error correctionprocessing of the present embodiment.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

(Configuration of Flow Velocimeter System)

A flow velocimeter system that is one embodiment of the presentinvention illustrated in FIG. 1 is to measure a moving mode of fluidflowing through a test section 50 (designated region). The flowvelocimeter system is used for designing and development of a movingbody such as a vehicle, to measure a moving mode of fluid around themoving body, for example. The flow velocimeter system includes alight-sheet generation device 10, an imaging device 20, an imageanalysis device 30 and a control unit 40.

The light-sheet generation device 10 includes a light source 11 and anoptical system 12 to form sheet-form light (light sheet) at the testsection 50. The light source 11 is configured to generate light(electromagnetic waves) such as laser light. The optical system 12 isconfigured so that light generated from the light source 11 is convertedinto sheet-form light, and then is applied at the test section 50. Theoptical system 12 includes a first cylindrical lens 12 a, a secondcylindrical lens 12 b and an output window 12 c.

The imaging device 20 has an optical axis that is disposed to beorthogonal to the light sheet (laser light sheet) formed at the testsection 50. The imaging device 20 forms an image of a group of particlesincluded in the fluid flowing through the test section 50intermittently. Each of a plurality of pixels making up the image formedhas a brightness resulting from the scattered light of the particles atthe light sheet, for example, as a pixel value.

The image analysis device 30 includes a computer (made up of anarithmetic processing unit or one or a plurality of processors (CPU), amemory (storage device), an I/O circuit and the like). The imageanalysis device 30 receives data of an image formed by the imagingdevice 20, and executes arithmetic processing for analysis of the image.The processor making up the image analysis device 30 reads necessarydata and such software from the memory, and executes arithmeticprocessing to the data in accordance with the software.

The control unit 40 includes a computer that is common to the computermaking up the image analysis device 30 or that is different from thecomputer. A generation timing of a light sheet by the light-sheetgeneration device 10 and an imaging timing by the imaging device 20 aresynchronized by various synchronization methods suitable for the type ofthe light source 11. As illustrated in FIG. 2, the control unit 40controls the light-sheet generation device 10 to intermittently generatea light sheet at the end of the exposure duration of the odd number oftimes by the imaging device 20 (2q+1-th, 2q+3-th, where q is an integer)and at the beginning of the exposure duration of the even number oftimes (2q+2-th, 2q+4-th). Such intermittent duration is adjusted so thatthe standard moving amount of a target at the formed image becomes abouta predetermined number of pixels (e.g., 5 pixels).

The test section 50 is defined with a pipe having a rectangular crosssection that is made of transparent acrylic, where the distance from thelens making up the imaging device 20 to the irradiation face of thelight sheet is set at Z₀. Fluid (air) including particles mixed thereinbeforehand is allowed to flow through the inside of the pipe. Light atthe light sheet is scattered by the particles, whereby an image ofparticles including pixels having a brightness resulting from thescattered light can be acquired.

An image of an independent particle at the image of particles is calleda “particle image”. An image of a group of particles made up ofcollected particles at the image of particles is called a “particlegroup image”. A three-dimensional orthogonal coordinate system(coordinate system defined with X-axis and Y-axis that are parallel tothe light sheet and Z-axis that is perpendicular to the light sheet)representing the location of a particle group at the test section 50 iscalled a “real-space coordinate system”. A two-dimensional orthogonalcoordinate system (coordinate system defined with x-axis and y-axis thatare parallel to the image-forming plane) representing the location of aparticle group, for example, in the particle image is called an “imagecoordinate system”.

Herein, the relationship between the coordinate (X, Y, Z) in thereal-space coordinate system and the coordinate (x, y) in the imagecoordinate system is represented by the following Expressions (1) to (3)using the magnification of the imaging M:

$\begin{matrix}{x = \frac{X}{M}} & (1) \\{y = \frac{Y}{M}} & (2) \\{{Z_{0} - \frac{\Delta\; Z}{2}} \leq Z \leq {Z_{0} + \frac{\Delta\; Z}{2}}} & (3)\end{matrix}$

The location of each pixel of the formed image is represented with apixel coordinate (i, j) having a one-to-one correspondence to eachpixel. Then, i, j is converted in to x, y as stated above based on thesize of each pixel, followed by conversion with Expressions (1) to (2),whereby the location in the real space corresponding to each pixel (i,j) can be found. Brightness of each pixel is represented as I(i, j).

(Functions of Flow Velocimeter System)

Referring to FIG. 3, the following describes the outline of functions ofthe flow velocimeter system.

The imaging device 20 forms an image of a particle group flowing throughthe test section 50 intermittently (FIG. 3/STEP010). The formed imageincludes high-brightness pixels resulting from the scattered light ofthe particle group at the light sheet that is formed at the test section50, and the existence of the particles is represented by such a part ofthe high-brightness pixels. Herein, the time interval of theintermittent imaging is adjusted to enable the observation of a movingmode of each particle group at the test section 50.

The image analysis device 30 executes noise reduction processing to eachformed image (FIG. 3/STEP020). The noise reduction processing includesambient light reduction processing (FIG. 3/STEP021) and particle imagecorrection processing (FIG. 3/STEP022).

The processing to reduce the influence from ambient light (FIG.3/STEP021) yields a subtraction image that is a result of thesubtraction of a pixel value of each pixel making up a background imagefrom a pixel value (brightness) of each pixel making up the formed imageat each time. The “background image” means an image at the test section50 that is formed by the imaging device 20 when a light sheet is formedthere but there are no particles. This can reduce the influence ofambient light (e.g., reflected light from a body of the vehicle) fromthe formed image at each time.

The particle image correction processing (FIG. 3/STEP022) corrects thebrightness of all pixels making up a predetermined range that is atleast a part of the formed image based on the brightness of a part ofpixels making up the predetermined range to correct the size of theimage of each particle group image. Specifically the minimum value ofthe brightness is subtracted from the brightness of pixels making up thepredetermined range.

FIG. 4 shows the relationship between the size of a particle image(particle image size) used when the PIV processing is executed to eachparticle image at different times that are created artificially usingrandom numbers and an error (average) for the correct value of thecorrelation value between inspection regions defined for the particleimages at the different times. The “inspection region” is defined withthe inspection window used at the local flow velocity vector acquisitionprocessing (FIG. 3/STEP030) described later. The horizontal axisrepresents the particle image size in the units of pixels of the image.The vertical axis represents the corresponding correlation value error(average). In three cases having difference widow lengths WL(s) of theinspection window (WL(s1)<WL(s2)<WL(s3)), the curves on the relationshipbetween the image size and the correlation value error are representedwith the broken line, the solid line and the dot-and-dash line.Typically a larger window length WL(s) of the inspection window means asmaller correlation value error between the inspection regions definewith such an inspection window.

As can be seen from FIG. 4, the correlation value error betweeninspection regions becomes the minimum when the particle image size isabout A pixel, irrespective of the window length WL(s). The “A pixel” is2-pixel plus for the cross correlation based on two images, for example.Presumably a too large particle image size makes it difficult todistinguish the particle image and the edge because the edge (borderregion) of one pixel image overlaps with the edge of another particleimage.

As illustrated in FIG. 5, the particle image correction processing (FIG.3/STEP022) subtracts the minimum brightness I_(min)(i₁, j₁) at the(A+1)×(A+1) pixels from the brightness I(i,j) of the pixels included inthe predetermined range having the size of (A+1)×(A+1) pixels (e.g., 3×3pixels). This results in the brightness at the edge of the particleimage to be zero (reference value), and so as illustrated in FIG. 5, theoriginal particle image size indicated with the broken line is correctedto be a particle image size (about A×A pixels) indicated with the solidline. This can minimize the correlation value error between inspectionregions that are defined in particle images at different times, andaccordingly leads to the improvement of measurement accuracy of the flowvelocity vector.

Specifically if the particle image correction processing is omitted, theoutline of each particle becomes blurred as illustrated in the upperpart of FIG. 6A. This makes the range in the formed image foracquisition of the flow velocity vector at the flow-velocity vectoracquisition processing (FIG. 3/STEP030) relatively narrower asillustrated in the lower part of FIG. 6A. On the other hand, when theparticle image correction processing is executed, then the outline ofeach particle image formed becomes clear as illustrated on the upperside of FIG. 6B. This means that the range in the formed image foracquisition of the flow velocity vector at the flow-velocity vectoracquisition processing (FIG. 3/STEP030) becomes relatively wider asillustrated on the lower side of FIG. 6B. Herein not the brightnessminimum value but the brightness average of all pixels included in thepredetermined range may be subtracted from the brightness of the pixelsincluded in the predetermined range.

Following the noise reduction processing (FIG. 3/STEP20), the imageanalysis device 30 performs processing to acquire the flow velocityvector v(a, b, T) from a plurality of images formed based on thecorrelation relationship between a first inspection region included in afirst image at a first time T and a second inspection region included ina second image at a second time T+ΔT (FIG. 3/STEP030). The time ΔT is asmall time (e.g., intermittent imaging time interval), and the firstimage and the second image are two images that are continuous in termsof clock time, for example.

Herein as illustrated in FIG. 7A, the brightness of the pixels includedin the first image is represented as I¹(i,j) using the superscriptnumber 1, the pixel coordinate value i of the pixel in x-axis directionand the pixel coordinate value j of the pixel in y-axis direction. Thenas illustrated in FIG. 7A, the brightness of the pixels included in thesecond image is represented as I²(i,j) using the superscript number 2and the pixel coordinates (i, j) of the pixel. When the types of theimage such as the first image and the second image are notdistinguished, it is simply represented as I(i,j).

In the first image, a plurality of first inspection regions R¹(a,b) aredefined. Herein, “a” denotes an index value indicating how they aredisplaced in the x-axis direction with reference to the inspectionregion located at the upper left corner of the first image. Then “b”denotes an index value indicating how they are displaced in the y-axisdirection with reference to the inspection region located at the upperleft corner of the first image. As illustrated in the range surroundedwith the bold line in FIG. 7A, any first inspection region R¹(a,b) is arectangle or a square of m pixels horizontally (x-axis direction) and npixels vertically (y-axis direction).

For each of the first inspection regions R¹(a,b), one first velocitydefining point p¹(a,b) is defined inside or on the border of the region.The first velocity defining point p¹(a,b) is located with reference tothe first inspection regions R¹(a, b) at a constant position (e.g., thecenter of the region) irrespective of the difference of the index value(a,b). The first velocity defining point p¹(a,b) becomes one of thereference points to calculate the flow velocity vector v(a, b, T). Onevelocity defining point and another velocity defining point adjacent tothe one velocity defining point are distant by d_(x) (or d_(y)) pixelsin the x-axis direction (or in the y-axis direction).

Next, referring to FIGS. 7A to 7D, the distribution of brightness ateach first inspection region is described below. The distribution ofbrightness at the first inspection region R¹(1,1) located at the leftupper corner of the first image can be represented by the followingExpression (4) using the brightness of pixels included in the firstinspection region R¹(1,1).I ¹(i,j) where 1≦i≦m,1≦j≦n  (4)

At the first inspection region R¹(1, 1), one first velocity definingpoint p¹(1,1) is defined.

For instance, as illustrated in FIG. 7B, a first velocity defining pointp¹(2,1) that is adjacent to the first velocity defining point p¹(1,1) inthe x-axis positive direction is distant from the first velocitydefining point p¹(1,1) by d_(x) pixels in the x-axis positive direction.That is, the distribution of brightness at the first inspection regioncorresponding to the first velocity defining point p¹(2,1) isrepresented by Expression (5).I ¹(i+d _(x) ,j) where 1≦i≦m,1≦j≦n  (5)

Then as illustrated in FIG. 7C, a first velocity defining point p¹(1,2)that is adjacent to the first velocity defining point p¹(1,1) in they-axis positive direction is distant from the first velocity definingpoint p¹(1,1) by d_(y) pixels in the y-axis positive direction. That is,the distribution of brightness at the first inspection regioncorresponding to the first velocity defining point p¹(1,2) isrepresented by Expression (6).I ¹(i,j+d _(y)) where 1≦i≦m,1≦j≦n  (6)

As illustrated in FIG. 7D, a first velocity defining point p¹(a₁,b₁) isdistant from the first velocity defining point p¹(1,1) by d_(x)(a₁−1)pixels in the x-axis positive direction and by d_(y)(b₁−1) pixels in they-axis positive direction. That is, the distribution of brightness atthe first inspection region corresponding to the first velocity definingpoint p¹(a₁,b_(i)) is represented by Expression (7).I ¹(i+d _(x)(a ₁−1),j+d _(y)(b ₁−1)) where 1≦i≦m,1≦j≦n  (7)

In the second image, a plurality of second inspection regionsR²(a,b,α,β) are defined corresponding to the plurality of firstinspection regions R¹(a,b) defined in the first image. Herein “a” and“b” are the same index values as “a” and “b” of the first inspectionregion R¹(a,b). As illustrated in FIG. 8A, “α” is an index valueindicating the distance between a first inspection region R¹(a,b) as atarget and a second inspection region R²(a,b,α,β) in the x-axisdirection (precisely the distance between the first velocity definingpoint p¹(a,b) and a second velocity defining point p²(a,b,α,β) describedlater in the x-axis direction) in the units of pixels. As illustrated inFIG. 8A, “β” is an index value indicating the distance between a firstinspection region R¹(a,b) as a target and a second inspection regionR²(a,b,α,β) in the y-axis direction (precisely the distance between thefirst velocity defining point p¹(a,b) and the second velocity definingpoint p²(a,b,α,β) described later in the y-axis direction) in the unitsof pixels. When the first inspection region R¹(a,b) and the secondinspection region R²(a,b) are not distinguished, they are represented asan inspection region R⁰(a,b) using R, the index values a, b and thesuperscript number 0.

Each of the second inspection regions R²(a,b,α,β) has the same shape andsize as those of the first inspection region R¹(a,b). That is, each ofthe second inspection regions R²(a,b,α,β) is a rectangular orsquare-shaped region of m pixels horizontally (x-axis direction) and npixels vertically (y-axis direction).

For each of the second inspection regions R²(a,b,α,β), one secondvelocity defining point p²(a,b,α,β) is defined inside or on the borderof the region. The second velocity defining point p²(a,b,α,β) is locatedwith reference to the second inspection regions R²(a, b) at a constantposition (e.g., the center of the region) irrespective of the differenceof the index value (a,b,α,β). The first inspection region R¹(a,b) andthe first velocity defining point p¹(a,b) have the same positionalrelationship as that of the second inspection region R²(a,b,α,β) and thesecond velocity defining point p²(a,b,α,β). Such a second velocitydefining point p²(a,b,α,β) becomes one of the reference points tocalculate the flow velocity vector v(a, b, T).

As indicated with the range surrounded with the bold line in FIG. 8A,the second inspection region R²(a,b,α,β) is a region of m pixelshorizontally (x-axis direction) and n pixels vertically (y-axisdirection) in size in the second image.

The distribution of brightness at the second inspection region R²(a,b,α,β) is represented by Expression (8).I ²(i+dx(α−1)+α,j+d(b−1)+β) where 1≦i≦m,1≦j≦n  (8)

For more detailed descriptions of the processing at STEP030, the flowvelocimeter system determines the size m, n of the inspection region bythe image analysis device 30 (FIG. 3/STEP031).

The flow velocimeter system increases m and n (e.g., m=64, n=64) at thetime of starting of the processing to perform correlation processing atSTEP032 in a wider range. Then as the processing proceeds, the flowvelocimeter system calculates the maximum movement amount of theparticles, x_max=max(|v(a, b, T_(x)|) in the x direction and the maximummovement amount of the particles, y_max=max(|v(a, b, T)_(y)|) in the ydirection, and sets the minimum m and n among m and n satisfyinginequalities (9) and (10). In other words, the flow velocimeter systemconsiders the maximum movement amount of the particles or the like andsets m and n so that the particles or the like after the movement areincluded in the searching region SR(a,b) (see Expression (13)) describedlater and the size of SR(a,b) becomes the minimum.

$\begin{matrix}{{- \frac{m}{2}} < {x\_ max} < \frac{m}{2}} & (9) \\{{- \frac{n}{2}} < {y\_ max} < \frac{n}{2}} & (10)\end{matrix}$

The subscript letters such as x and y represent the values of thecomponents as indicated in Expressions (11) and (12).v(a,b,t)_(x): value of x component of local flow velocity vector v  (11)v(a,b,t)_(xy): value of x component or y component of local flowvelocity vector v  (12)

As a result, as the processing proceeds, unnecessary calculation can bereduced because the size of the inspection region approaches the optimumsize in view of the movement amount of the particles or the like.

The image analysis device 30 calculates a correlation function C(α,β) ofthe brightness distribution of the first inspection region in the firstimage, I¹(i+d(a−1),j+d(b−1)) (1≦i≦m, 1≦j≦n) and the brightnessdistribution of the second inspection region in the second image,I²(i+d(a−1)+α, j+d(b−1)+β) ((1≦i≦m, 1≦j≦n) for all of the firstinspection regions R¹(a,b) included in the first image in the range of−m/2<α<m/2 and −n/2<β<n/2 (FIG. 3/STEP032).

More specifically, as illustrated in FIG. 8B, the flow velocimetersystem finds the correlation function C(α,β) having, as the output,correlation values c(α₁, β₁), c(α₂, β₂) and c(α₃, β₃) . . . between acertain first inspection region R¹(a₂,b₂) and the second inspectionregions R²(a₂, b₂, α₁, β₁), R²(a₂, b₂, α₂, β₂) and R²(a₂, b₂, α₃, β₃)included in the searching region SR(a₂,b₂) that is determinedcorresponding to the first inspection region R¹(a₂,b₂).

Herein, the searching region SR(a₂,b₂) is the range that can berepresented by Expression (13), for example.

$\begin{matrix}{{{{SR}\left( {a_{2},b_{2}} \right)} = {{{\left( {{i + {{dx}\left( {a_{2} - 1} \right)}},{j + {{dy}\left( {b_{2} - 1} \right)}}} \right)\mspace{14mu}{where}} - \frac{m}{2}} \leq i \leq \frac{3m}{2}}},{{- \frac{n}{2}} \leq j \leq \frac{3n}{2}}} & (13)\end{matrix}$

Herein, the correlation function C(α,β) is represented by Expression(14), for example.

$\begin{matrix}{{C\left( {\alpha,\beta} \right)} = {\sum\limits_{i = 1}^{m}{\sum\limits_{j = 1}^{n}{{I^{1}\left( {{i + {d\left( {a_{2} - 1} \right)}},{j + {d\left( {b_{2} - 1} \right)}}} \right)}{I^{2}\left( {{i + {d\left( {a_{2} - 1} \right)} + \alpha},{j + {d\left( {b_{2} - 1} \right)} + \beta}} \right)}}}}} & (14)\end{matrix}$

For the calculation of the correlation function C(α,β), FFT processingis performed to the first inspection region R¹(a₂,b₂) and the secondinspection region R²(a₂,b₂,α,β) for each of the vertical direction (ydirection) and the horizontal direction (x direction) as illustrated inFIG. 9. Then, the brightness distribution characteristic of the firstinspection region R¹(a₂,b₂) and the brightness distributioncharacteristic of the second inspection region R²(a₂,b₂,α,β), which areacquired by the FFT processing, are compared. In this way, thecorrelation function C(α,β) can be calculated.

Since the adjacent inspection regions are overlapped partially, the FFTprocessing result for such an overlapping part is preferably reused. Forinstance, as illustrated in FIG. 10A, the flow velocimeter system reusesthe FFT processing result at the inspection region R⁰(a₃,b₃) orR⁰(a₄,b₄) for the FFT processing of the inspection region R⁰(a₃+1,b₃) orR⁰(a₄,b₄+1)) adjacent to the inspection region R⁰(a₃,b₃) or R⁰(a₄,b₄).

As illustrated in FIG. 10B, after the FFT processing is executed for theinspection region R⁰(a₃,b₃), FFT processing is executed as follows forthe inspection region R⁰(a₃+1,b₃) that is displaced to the right (+xdirection) by 0.25 times the width of the inspection region R⁰(a₃,b₃).The FFT processing in the y direction (see left side of FIG. 10B) cancalculate the y-direction component of the flow velocity vector and theFFT processing in the x direction (see right side of FIG. 10B) cancalculate the x-direction component of the flow velocity vector. In thiscase, the result of the FFT processing in the y direction relating tothe overlapping region (accounting for 75% of the entire) with theinspection region R⁰(a₃,b₃) is applied to the inspection regionR⁰(a₃+1,b₃), whereby the FFT processing in the y direction of such anoverlapping region can be omitted.

As illustrated in FIG. 10C, after the FFT processing is executed for theinspection region R⁰(a₄,b₄), FFT processing is executed as follows forthe inspection region R⁰(a₄+1,b₄) that is displaced downward (+ydirection) by 0.25 times the height of the inspection region R⁰(a₄,b₄).The FFT processing in the y direction (see left side of FIG. 10C) cancalculate the y-direction component of the flow velocity vector and theFFT processing in the x direction (see right side of FIG. 10C) cancalculate the x-direction component of the flow velocity vector. In thiscase, the result of the FFT processing in the x direction relating tothe overlapping region (accounting for 75% of the entire) with theinspection region R⁰(a₄,b₄) is applied to the inspection regionR⁰(a₄,b₄+1), whereby the FFT processing in the x direction of such anoverlapping region can be omitted.

Such reusing of the FFT processing result can reduce the amount ofcalculation during the processing roughly represented by Expression (15)as compared with the case of not reusing the FFT processing result.

$\begin{matrix}{\frac{1}{2}\left( {\frac{m - d_{x}}{m} + \frac{n - d_{y}}{n}} \right)} & (15)\end{matrix}$

The flow velocimeter system finds the set of α₀ and β₀ such that thethus found correlation function C(α,β) yields the maximumc_(max)(α₀,β₀). The second inspection region R²(a₂,b₂,α₀,β₀) yieldingthe maximum correlation function C(α,β) with the first inspection regionR¹(a₂, b₂) can be considered to include the image of the particles thatis the same or has the same arrangement as those of the image of theparticles included in the first inspection region R¹(a₂, b₂).

When α₀=a₃, β₀=β₃, as illustrated with the moving vector D(a,b) in FIG.8C, the particles or the like existing on the real space coordinatesystem corresponding to the first inspection region R¹(a₂,b₂) at thetime T presumably moves to the position on the real space coordinatesystem corresponding to the second inspection region R²(a₂, b₂, α₃, β₃)by the time T+ΔT.

That is, the flow velocity vector v(a₂, b₂, T) in the first inspectionregion R¹(a₂,b₂) at the time T can be found by Expression (16).

$\begin{matrix}{{v\left( {a_{2},b_{2},T} \right)} = {\left( \frac{D\left( {a_{2},b_{2}} \right)}{\Delta\; T} \right) = \begin{pmatrix}\frac{\alpha_{0}}{\Delta\; T} \\\frac{\beta_{0}}{\Delta\; T}\end{pmatrix}}} & (16)\end{matrix}$

Repeating the procedure at STEP032 allows the flow velocity vector v(a,b, T) at the time T in each of the first inspection regions R¹(a,b) tobe acquired.

The flow velocimeter system makes the image analysis device 30 toperform the processing to correct the flow velocity vector v(a, b, T)that can be considered as a measurement error, while considering theflow velocity vector v(a_(±), b_(±), T) of the first inspection regionR¹(a_(±),b_(±)) around the first inspection region R¹(a,b) on the basisof the flow velocity vector v(a, b, T) in the first inspection regionR¹(a,b) at the time T (FIG. 3/STEP040).

More specifically, determination is made as to whether each componentv(a, b, T)_(xy) of the flow velocity vector v(a, b, T) in the firstinspection region R¹(a,b) at the time T satisfies a condition (a spatialerror correction condition) that is determined from each component v(a,b, T)_(xy) of the flow velocity vector v(a_(±), b_(±), T) of thesurrounding first inspection region R¹(a_(±), b_(±)). The spatial errorcorrection condition that can be used includes that each component v(a,b, T)_(xy) of the flow velocity vector v(a, b, T) is not included in apredetermined range (e.g., the range that is determined usingdispersion) that is determined from the median value or the average ofeach component v(a_(±), b_(±), T)_(xy) of the flow velocity vectorv(a_(±), b_(±), T) in the surrounding first inspection region R¹(a_(±),b_(±)).

When the flow velocity vector v(a, b, T) satisfies the spatial errorcorrection condition, each component v(a, b, T)_(xy) of the flowvelocity vector of the first inspection region R¹(a_(±), b_(±)) isreplaced with each component v(a_(±), b_(±), T)_(xy) of the surroundingflow velocity vector. For instance, the weighted average of thecomponents v(a_(±), b_(±), T)_(xy) of the surrounding flow velocityvector is used as each component v(a, b, T)_(xy) of the flow velocityvector of the first inspection region R¹(a_(±), b_(±)).

The flow velocimeter system determines whether the image analysis device30 performs the processing of STEP030 to STEP040 to all images or not(or whether there are images at the next time (third time and fourthtime) of the first time T and the second time T+ΔT or not) (FIG.3/STEP050). When the determination result is negative (FIG. 3/STEP050 .. . NO), the flow velocimeter system continues the processing of STEP030to STEP040 while setting the third time as a new first time T.

The processing at FIG. 3/STEP030 to STEP050 yields the flow velocityvector v(a, b, T) in the first inspection region R¹(a,b) at the time Tas indicated with the solid line in FIG. 12A.

When the determination result is affirmative (FIG. 3/STEP050 . . . YES),the flow velocimeter system performs time-series error correctionprocessing by the image analysis device 30 (FIG. 3/STEP060).

The time-series error correction processing is to correct a target flowvelocity vector v(a, b, T) in the first inspection region R¹(a,b) at thetime T indicated with the solid line in FIG. 12A to a target flowvelocity vector v₊(a, b, T) in the first inspection region R¹(a,b) atthe time T indicated with the dashed line in FIG. 12A while consideringa reference flow velocity vector v(a, b, T_(±)) in the first inspectionregion R¹(a,b) at the time T_(±) before and after the time T.

Herein, T−WL(s)/2≦T_(±)≦T+WL(s)/2, and WL(s) denotes the window length.

Following that, the flow velocimeter system performs the processing forthe component v(a, b, T)_(xy) of each axis (e.g., x axis of the flowvelocity vector of the particle group) of the flow velocity vector v(a,b, T).

The flow velocimeter system initializes a certain time T (FIG.11/STEP061).

The flow velocimeter system set the window at the certain time T (FIG.11/STEP062).

When the window length WL(s) is 5, then as illustrated in FIG. 13, thefive flow vectors (in the case of T=t₃, v(a, b, t₁) to v(a, b, t₅))before and after the time T (in the case of T=t₃, time t₁, t₂, t₃, t₄and t₅) are used in the following processing.

Next, the flow velocimeter system defines various types of parametersfor condition determination (FIG. 11/STEP063). Specifically theparameters are defined by Expressions (17) and (18).v_rms_(xy): root mean square (rms) of each component of local flowvelocity vector in window  (17)v_mid_(xy): median value of each component of each local flow velocityvector in window  (18)

Next, the flow velocimeter system determines whether the target flowvelocity vector v(a, b, T) in the first inspection region R¹(a,b) attime T satisfies the time-series error correction condition or not (FIG.11/STEP064). Specifically the determination is made as to whetherinequality (19) is satisfied or not.(v(a,b,t)_(xy) −v_mid_(xy))²>2v_rms_(xy) ²  (19)

Herein, the left side of the inequality (19) corresponds to the “degreeof difference” of the present invention, and the right side of theinequality (19) corresponds to the “threshold” of the present invention.In this way, the threshold is determined corresponding to the value ofcomponent of each flow velocity vector in the window, whereby the degreeof difference of the flow velocity vector can be determined flexibly onthe basis of the mode of the flow velocity vector at the times beforeand after the time t as compared with the case where the threshold isdetermined beforehand.

The degree of difference and the threshold are not limited to the aboveexpressions, and a difference between v(a, b, T)_(xy) and the average ofother flow velocity vectors in the window, for example, may be used asthe degree of difference, and a value that is determined by anexperiment beforehand may be used as the threshold.

When the result of the determination is No (FIG. 11/STEP064 . . . NO),the flow velocimeter system does not correct the target flow velocityvector v(a, b, T)(velocity data) at the time T and proceeds to the nextprocessing.

For instance, as illustrated in FIG. 12A and FIG. 13, when it isdetermined that the target flow velocity vector v(a, b, t₃) at the timet₃ does not satisfy the inequality (19) while considering the referenceflow velocity vectors v(a, b, t₁), v(a, b, t₂), v(a, b, t₄) and v(a, b,t₅) at t₁, t₂, t₄ and t₅ before and after the time t₃, the flowvelocimeter system does not change the value of the flow velocity vectorv(a, b, t₃) and proceeds to the processing at the next time t₄.

When the result of the determination is Yes (FIG. 11/STEP064 . . . YES),the flow velocimeter system corrects the target flow velocity vectorv(a, b, T) using the weighted average of each component of the referenceflow velocity vectors v(a, b, T_(±)) at other times in the window (e.g.,by increasing the weighting coefficient of the reference flow velocityvectors at times T−1 and T+1 close to the time T and decreasing theweighting coefficient of the reference flow velocity vectors at timesT−2 and T+2 that are away from the time T) (FIG. 11/STEP065).Specifically the target flow velocity vector v(a, b, T) is corrected byExpression (20), for example.

$\begin{matrix}\left. {v\left( {a,b,T} \right)}_{xy}\leftarrow\frac{\begin{matrix}{{\frac{1}{2}{v\left( {a,b,{T - 2}} \right)}_{xy}} + {v\left( {a,b,{T - 1}} \right)}_{xy} +} \\{{v\left( {a,b,{T + 1}} \right)}_{xy} + {\frac{1}{2}{v\left( {a,b,{T + 2}} \right)}_{xy}}}\end{matrix}}{\frac{1}{2} + 1 + 1 + \frac{1}{2}} \right. & (20)\end{matrix}$

For instance, when it is determined that the target flow velocity vectorv(a, b, t₄) at the time t₄ satisfies the inequality (19) whileconsidering the reference flow velocity vectors v(a, b, t₂), v(a, b,t₃), v(a, b, t₅) and v(a, b, t₆) at t₂, t₃, t₅ and t₆ before and afterthe time t₄, the flow velocimeter system corrects the flow velocityvector v(a, b, t₄) to v₊(a, b, t₄) by Expression (20) as indicated withthe arrow Q1 in FIG. 12A.

Then, the flow velocimeter system uses the flow velocity vector v₊(a, b,t₄) after correction in the processing at the time t₄ or later. Forinstance, when it is determined that the flow velocity vector v(a, b,t₅) at the time t₅ satisfies the inequality (19) while considering thereference flow velocity vectors v(a, b, t₃), v(a, b, t₄) (aftercorrection), v(a, b, t₆) and v(a, b, t₇) at t₃, t₄, t₆ and t₇ before andafter the time t₅, the flow velocimeter system corrects the flowvelocity vector v(a, b, t₅) to v₊(a, b, t₅) by Expression (20) asindicated with the arrow Q2 in FIG. 12A.

The flow velocimeter system may determine before the processing at FIG.11/STEP065 whether the ratio of the reference flow velocity vectors v(a,b, T_(±)) satisfying the condition represented by the inequality (19) tothe entire reference flow velocity vectors v(a, b, T_(±)) (hereinafterthis is called “satisfying ratio of inequality (19)”) is a predeterminedratio (e.g., 40%) or more. When the satisfying ratio of inequality (19)is the predetermined ratio or more, the flow velocimeter systemincreases the window length WL(s) until the satisfying ratio ofinequality (19) becomes less than the predetermined ratio, and then mayperform the processing at FIG. 11/STEP063 to STEP065. As illustrated inFIG. 12B, this can increase the ratio of the reference flow velocityvectors that are likely to be accurate due to the increased windowlength WL(s) when the satisfying ratio of inequality (19) of thereference flow velocity vectors is high, i.e., when there are manyreference flow velocity vectors that are likely to be noise. As aresult, the target flow velocity vector can be corrected more suitablyfrom the viewpoint of bringing it to the correct value.

The flow velocimeter system determines whether the processing isperformed for all of the times T or not (FIG. 11/STEP066). When theresult of the determination is No, the flow velocimeter system increasesT by 1 (FIG. 11/STEP067) and repeats the processing at STEP062 toSTEP065. When the result of the determination is Yes, the flowvelocimeter system finishes the time-series error correction processing.

The flow velocimeter system may be configured so as to determine, beforeor after the processing at FIG. 11/STEP066 at time T, whether thecorrection ratio of flow velocity vectors in a part or all of the pastduration (e.g., time T−γ to time T(γ>0)) including the time T is apredetermined value or more (e.g., the correction ratio of the flowvelocity vector during the duration is 5% or more), and when the resultof the determination is Yes, to execute the processing at STEP061 toSTEP067 during the duration repeatedly. When the correction ratio of theflow velocity vectors is high, original data includes a lot of flowvelocity vectors (noise) being likely to be a measurement error, and sothe data after correction also is likely to include a lot of noiseremaining. In such a case, the processing at STEP061 to STEP067repeatedly executed can decrease the remaining ratio of the noise.

Even when the correction ratio of flow velocity vectors during certainduration is a predetermined value or less, this may be the case wherethe correction ratio of flow velocity vectors during duration shorterthan the certain duration is high, i.e., where noise concentrateslocally.

In view of this, the flow velocimeter system may be configured so thatthe predetermined value increases with decreases in the length of theduration (e.g., the predetermined value is a first predetermined value(e.g., 5%) during past first duration (time T−γ1 to time T(γ1>0)), andis a second predetermined value (e.g., 10%) higher than the firstpredetermined value during second duration (time T−γ2 to timeT(γ1>γ2>0)) shorter than the first duration. The first duration may bethe past entire duration.

As a result, the processing at STEP061 to STEP067 is executed repeatedlyonly during the duration having an especially high correction ratio offlow velocity vectors, and so the probability of local concentration ofnoise can be decreased while suppressing an increase in the load for theprocessing.

After the time-series error correction processing (FIG. 3/STEP060), theflow velocimeter system finishes a series of the processing.

(Advantageous Effects of the Flow Velocimeter System)

According to the thus configured flow velocimeter system, a light sheetis generated at a designated region, and images of fluid flowing throughthe designated region are formed at a plurality of different times. Foreach of a plurality of inspection regions making up the images, when aflow velocity vector v(a, b, T) at a certain time T is different from aflow velocity vector (reference flow velocity vectors) v(a, b, T_(±)) atone or a plurality of times T_(±) that is different from the certaintime T beyond a threshold v_rms, the flow velocity vector v(a, b, T) atthe certain time T is corrected (FIG. 11/STEP064 . . . YES→STEP065, seeinequality (19) and Expression (20)). As a result, a flow velocityvector v(a, b, T) that is likely to be noise different from the originalform can be corrected closer to the original form. Then, the accuracy ofmeasurement of time-series or changing mode of the fluid velocity can beimproved.

Specifically FIG. 14A illustrates original data (dotted line) and databefore the time-series error correction processing (solid line). FIG.14B illustrates original data (dotted line) and data after thetime-series error correction processing by the technique described inPatent Document 1 (solid line). FIG. 14C illustrates original data(dotted line) and data after the time-series error correction processingaccording to the present embodiment (solid line).

As FIG. 14A to FIG. 14C show, the data after the time-series errorcorrection processing by the conventional technique (see FIG. 14B) isdecreased in a part of the data that is extremely off from the originaldata, but the data as a whole is rather different from the originaldata. On the other hand, the data after the time-series error correctionprocessing according to the present embodiment (solid line) includeserror data decreased, and the reproducibility of the original data canbe increased as a whole. That is, the flow velocimeter system that isone embodiment of the present invention can improve the time-seriesmeasurement accuracy of the fluid velocity.

Brightness of all pixels making up a predetermined range that is atleast a part of the formed image is corrected (FIG. 3/STEP022, see FIG.4 and FIG. 5). This enables clear discrimination of the brightnessresulting from the border between particles included in the fluid ormaking up the fluid and the background light and the brightnessresulting from the scattered light of the particles (see FIG. 6B upperside). That is, the measurement accuracy of a time-series correlation ofthe region including the particles in the formed image can be improved,and accordingly the time-series measurement accuracy of the fluidvelocity can be improved (see FIG. 6B lower side).

For the calculation of the correlation function C(α,β), when there is anoverlapped region in one inspection region and in another inspectionregion among the plurality of inspection regions, a result of the fastFourier transformation processing at the other inspection region isapplied to the overlapped region for the fast Fourier transformation atthe one inspection region (see FIG. 9A and FIG. 9B as well as FIG. 10Ato FIG. 10C). This can save the calculation resource for the fastFourier transformation processing (FFT processing) in the overlappedregion between the plurality of inspection regions. This means that thefluid velocity having a high pixel density of the image also can bemeasured effectively in a time-series manner, and so the measurementaccuracy can be further improved.

The threshold v_rms can be determined adaptively corresponding to thereference flow velocity vector v(a, b, T_(±)), and so the time-seriesmeasurement accuracy of the fluid velocity can be further improved.

When the occupancy of different flow velocity vectors (reference flowvelocity vectors that do not satisfy the inequality (19)) among thereference flow velocity vectors v(a, b, T_(±)) at the plurality of timesexceeds a predetermined value, the number of the reference flow velocityvectors v(a, b, T_(±)) are corrected for increasing by correcting thewindow length WL(s) for increasing so that the occupancy becomes thepredetermined value or less. Then, the flow velocity vector v(a, b, T)at each time T can be corrected suitably based on the reference flowvelocity vector v(a, b, T_(±)) that is likely not to correspond to noiseor to correspond to the true value, and so the time-series measurementaccuracy of the fluid velocity can be improved.

What is claimed is:
 1. A flow velocimeter system comprising: alight-sheet generation device; an imaging device; and an image analysisdevice including a processor, the image analysis device being configuredto measure, based on images of fluid flowing through a designated regionthat are formed a plurality of times at different times by the imagingdevice while the light-sheet generation device generates a light sheetcontinuously or intermittently at the designated region, a time-serieschanging mode of a local flow velocity vector of the fluid at each of aplurality of inspection regions defined at the images in accordance witha cross correlation method, wherein the image analysis device isconfigured to, for an inspection region of the plurality of inspectionregions that has a degree of difference exceeding a threshold betweenthe local flow velocity vector at a certain time and a reference flowvelocity vector as the local flow velocity vectors at one or a pluralityof times that is different from the certain time, correct the local flowvelocity vector at the certain time with the reference flow velocityvector, wherein the image analysis device is configured to, whenoccupancy of different flow velocity vectors among the reference flowvelocity vectors at the plurality of times exceeds a predeterminedvalue, increase the number of the reference flow velocity vectors sothat the occupancy becomes the predetermined value or less, and whereinthe different flow velocity vectors are the reference flow velocityvectors at one time of the plurality of times, which have a degree ofdifference from the reference flow velocity vectors at another time thatexceeds the threshold.
 2. The flow velocimeter system according to claim1, wherein the image analysis device is configured to correct brightnessof all pixels making up a predetermined range that is at least a part ofthe images based on brightness of a pixel as a part making up thepredetermined range, and then calculate the local flow velocity vectorat each of the plurality of inspection regions.
 3. The flow velocimetersystem according to claim 1, wherein the image analysis device isconfigured to execute fast Fourier transformation processing at each ofthe plurality of inspection regions, thus calculating the local flowvelocity vectors, and when there is an overlapped region in oneinspection region and in another inspection region among the pluralityof inspection regions, to apply a result of the fast Fouriertransformation processing at the other inspection region to theoverlapped region for the fast Fourier transformation at the oneinspection region.
 4. The flow velocimeter system according to claim 1,wherein the image analysis device uses a value that is determined inaccordance with the reference flow velocity vector as the threshold.